pacman::p_load(
  rio, # import funcs
  sf, # work with spatial data
  here, # create relative paths
  janitor, # data cleaning
  lubridate, # date handling
  tidyverse # data science
)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")

library(tidyverse)

# linelist

path_ll <- fs::dir_ls(here::here("data", "final"), 
                      regex = "msf_linelist") |> max()

dat_raw <- import(path_ll) |> as.tibble()

# lab data
path_lab <- fs::dir_ls(here::here("data", "final"), 
                       regex = "msf_laboratory") |> max()

lab_raw <- import(path_lab) |> as_tibble()

# admin data
admin_1 <- st_read(here::here("data", "gpkg", "GEO-EXPORT-TCD-2024-04-11.gpkg"), layer = "ADM1", quiet= TRUE)
admin_2 <- st_read(here::here("data", "gpkg", "GEO-EXPORT-TCD-2024-04-11.gpkg"), layer = "ADM2", quiet = TRUE)

Introduction

This document provides directions for the analysis of the fake Measles dataset msf_linelist_moissala_2023-09-24.xlsx and the its corresponding laboratory dataset msf_laboratory_moissala_2023-09-24.xlsx.

Data Cleaning

Age classification

For measles outbreaks, it makes sense to use the following age classification:

0 - 11 months
< 6 months
6 - 8 months
9 - 11 months
1 - 4 years
5 - 14 years
15+ years

Middle Upper Arm Circumferences (MUAC)

MUAC is classified as follow:

Green (125+ mm)
Yellow (115 - 124 mm)
Red (<115 mm)

Epidemiological classification

confirmed: cases with a positive PCR result in lab data probable: cases with fever, coughand a rash suspected: all other cases

dat <- dat_raw |>
  # standardise variable names
  janitor::clean_names() |>
  # manually rename
  rename(
    id = epi_id_number,
    sex = sex_patient,
    age_unit = age_units_months_years,
    date_onset = date_of_onset_of_symptoms,
    hospitalisation = hospitalisation_yes_no,
    date_admission = date_of_admission_in_structure,
    date_outcome = date_of_outcome,
    sub_prefecture = sub_prefecture_of_residence,
    region = region_of_residence,
    fever = participant_had_fever,
    rash = participant_had_rash,
    cough = participant_had_cough,
    red_eye = participant_had_red_eye,
    pneumonia = participant_had_pneumonia,
    encephalitis = participant_had_encephalitis,
    muac = middle_upper_arm_circumference_muac,
    vacc_status = vaccination_status,
    vacc_doses = vaccination_dosage,
    outcome = patient_outcome_death_recovered_lama,
    site = msf_site,
    malaria_rdt = malaria_rdt
  ) |>
  # recoding
  mutate(
    sex = case_when(
      sex %in% c("f", "femme") ~ "female",
      sex %in% c("m", "homme") ~ "male",
      .default = sex
    ), 
    across(contains("date_"), ~ as.Date(.x)),
    across(
      c(
        fever,
        rash,
        cough,
        red_eye,
        pneumonia,
        encephalitis
      ),
      ~ case_match(.x, 
                   "Yes" ~ TRUE, 
                   "No" ~ FALSE, .default = NA)
    )
  ) |>
  # Categorise variables
  mutate(
    age_group = case_when(
      age_unit == "months" & age < 6 ~ "< 6 months",
      age_unit == "months" & between(age, 6, 8) ~ "6 - 8 months",
      age_unit == "months" & between(age, 9, 11) ~ "9 - 11 months",
      age_unit == "years" & between(age, 1, 4) ~ "1 - 4 years",
      age_unit == "years" & between(age, 5, 14) ~ "5 - 14 years",
      age_unit == "years" & between(age, 15, 40) ~ "15+ years"
    ),
    age_group = fct_relevel(
      age_group,
      c(
        "< 6 months",
        "6 - 8 months",
        "9 - 11 months",
        "1 - 4 years",
        "5 - 14 years",
        "15+ years"
      )
    ),
    muac_cat = case_when(
      muac >= 125 ~ "Green (125+ mm)",
      between(muac, 115, 124) ~ "Yellow (115 - 124 mm)",
      muac < 115 ~ "Red (<115 mm)"
    )
  ) |>
  relocate(
    age_group,
    .after = age_unit
  ) |>
  relocate(muac_cat, .after = muac) |> 
  
  filter(date_onset >= as.Date("2023-01-01"))

Laboratory data

lab_clean <- lab_raw |>
  clean_names() |>
  rename(
    case_id = msf_number_id,
    lab_id = laboratory_id,
    date_test = date_of_the_test,
    test_result = final_test_result
  ) |>
  mutate(
    date_test = ymd(date_test),
    ct_value = round(ct_value, digits = 1)
  )

There are some duplicates in the laboratory results. Some case_id were tested multiple times if there was a inconclusive test_result. We need to find them, and take the last sample tested

reactable::reactable(lab_clean |> get_dupes(case_id))

lab_clean <- lab_clean |>
  filter(
    .by = case_id,
    date_test == max(date_test, na.rm = TRUE)
  )

Some samples were negative, so these cases are not cases and need to be removed from analysis

lab_clean |> count(test_result)

## # A tibble: 0 × 2
## # ℹ 2 variables: test_result <chr>, n <int>

We join the lab_clean to the main linelists using the case_id as key. Then remove the negative case, and create an epidemiological classification

dat <- left_join(dat, lab_clean, by = c("id" = "case_id"))

dat <- dat |>
  filter(is.na(test_result) | test_result == "positive") |>
  mutate(
    epi_cat = case_when(
      test_result == "positive" ~ "confirmed",
      rash == TRUE & fever == TRUE & cough == TRUE ~ "probable",
      .default = "suspected"
    ),
    epi_cat = fct_relevel(epi_cat, c("confirmed", "probable", "suspected"))
  )

reactable::reactable(dat |> tabyl(epi_cat) |> mutate(percent = round(percent * 100, digits = 2)))

Person

Demographics

dat |>
  select(
    sex,
    age_group,
    muac_cat,
    vacc_status
  ) |>
  gtsummary::tbl_summary(label = list(
    sex ~ "Gender",
    age_group ~ "Age groups",
    muac_cat ~ "MUAC category",
    vacc_status = "Vaccination status",
    malaria_rdt = "Malaria RDT",
    outcome = "Outcome"
  ))

Characteristic	N = 4,510¹
Gender
female	2,293 (51%)
male	2,217 (49%)
Age groups
< 6 months	1,188 (27%)
6 - 8 months	572 (13%)
9 - 11 months	441 (10%)
1 - 4 years	1,085 (25%)
5 - 14 years	847 (19%)
15+ years	263 (6.0%)
Unknown	114
MUAC category
Green (125+ mm)	3,739 (83%)
Red (<115 mm)	269 (6.0%)
Yellow (115 - 124 mm)	502 (11%)
Vaccination status
No	2,560 (68%)
Uncertain	520 (14%)
Yes - card	39 (1.0%)
Yes - oral	636 (17%)
Unknown	755
¹ n (%)

By sites

dat |>
  select(
    sex,
    age_group,
    muac_cat,
    vacc_status,
    site
  ) |>
  gtsummary::tbl_summary(
    by = site,
    label = list(
      sex ~ "Gender",
      age_group ~ "Age groups",
      muac_cat ~ "MUAC category",
      vacc_status = "Vaccination status",
      malaria_rdt = "Malaria RDT",
      outcome = "Outcome"
    )
  )

Characteristic	Bedaya Hospital, N = 839¹	Bekourou Hospital, N = 101¹	Bouna Hospital, N = 715¹	Danamadji Hospital, N = 57¹	Koumogo Hospital, N = 3¹	Moïssala Hospital, N = 2,795¹
Gender
female	418 (50%)	54 (53%)	360 (50%)	28 (49%)	2 (67%)	1,431 (51%)
male	421 (50%)	47 (47%)	355 (50%)	29 (51%)	1 (33%)	1,364 (49%)
Age groups
< 6 months	210 (26%)	32 (32%)	202 (29%)	12 (21%)	0 (0%)	732 (27%)
6 - 8 months	109 (13%)	16 (16%)	96 (14%)	4 (7.1%)	1 (33%)	346 (13%)
9 - 11 months	80 (9.7%)	11 (11%)	57 (8.2%)	6 (11%)	0 (0%)	287 (11%)
1 - 4 years	186 (23%)	23 (23%)	172 (25%)	19 (34%)	1 (33%)	684 (25%)
5 - 14 years	186 (23%)	17 (17%)	129 (19%)	12 (21%)	1 (33%)	502 (18%)
15+ years	50 (6.1%)	2 (2.0%)	38 (5.5%)	3 (5.4%)	0 (0%)	170 (6.2%)
Unknown	18	0	21	1	0	74
MUAC category
Green (125+ mm)	697 (83%)	85 (84%)	581 (81%)	53 (93%)	3 (100%)	2,320 (83%)
Red (<115 mm)	49 (5.8%)	4 (4.0%)	55 (7.7%)	2 (3.5%)	0 (0%)	159 (5.7%)
Yellow (115 - 124 mm)	93 (11%)	12 (12%)	79 (11%)	2 (3.5%)	0 (0%)	316 (11%)
Vaccination status
No	474 (68%)	56 (65%)	408 (69%)	27 (68%)	1 (50%)	1,594 (68%)
Uncertain	113 (16%)	13 (15%)	72 (12%)	6 (15%)	1 (50%)	315 (13%)
Yes - card	7 (1.0%)	0 (0%)	6 (1.0%)	1 (2.5%)	0 (0%)	25 (1.1%)
Yes - oral	107 (15%)	17 (20%)	104 (18%)	6 (15%)	0 (0%)	402 (17%)
Unknown	138	15	125	17	1	459
¹ n (%)

Age Pyramids

dat |>
  select(
    sex,
    age_group,
    site
  ) |>
  apyramid::age_pyramid(
    age_group = "age_group",
    split_by = "sex",
    proportional = TRUE,
    show_midpoint = TRUE
  ) +
  
  theme_minimal()

CFR analysis by site

# CFR only on known outcomes
dat |>
  summarise(
    .by = site,
    n_cases = n(),
    n_confirmed = sum(epi_cat == "confirmed"),
    n_deaths = sum(outcome == "dead", na.rm = TRUE),
    cfr = round(digits = 2, n_deaths / sum(outcome %in% c("recovered", "dead")) * 100)
  ) |>
  reactable::reactable(columns = list(
    n_cases = reactable::colDef(name = "N cases"),
    n_confirmed = reactable::colDef(name = "N confirmed"),
    n_deaths = reactable::colDef(name = "N deaths"),
    cfr = reactable::colDef(name = "CFR (%)")
  ))

Risks Factor analysis

Investigating age_group, muac_cat and vacc_status as risks factors for death

# Prepare the data for fitting the logistic regression
prep_logit <- dat |>
  # change group order for references
  mutate(
    age_group = fct_relevel(
      age_group,
      c(
        "15+ years",
        "< 6 months",
        "6 - 8 months",
        "9 - 11 months",
        "1 - 4 years",
        "5 - 14 years"
      )
    ),
    muac_cat = fct_relevel(
      muac_cat,
      c(
        "Green (125+ mm)",
        "Yellow (115 - 124 mm)",
        "Red (<115 mm)"
      )
    ),
    vacc_status = case_match(
      vacc_status,
      "Yes - card" ~ "Yes",
      "Yes - oral" ~ "Yes",
      "Uncertain" ~ NA,
      .default = vacc_status
    ),
    vacc_status = fct_relevel(
      vacc_status,
      c(
        "No",
        "Yes"
      )
    ),
    
    # outcome needs to be 1/0
    outcome_binary = case_when(
      outcome == "recovered" ~ 0,
      outcome == "dead" ~ 1,
      .default = NA
    )
  )

# fit the logistic regression
mdl <- glm(outcome_binary ~ sex + age_group + vacc_status + muac_cat, data = prep_logit, family = "binomial")

# view coeff
gtsummary::tbl_regression(
  mdl,
  exp = TRUE,
  label = list(
    sex ~ "Gender",
    age_group ~ "Age groups",
    muac_cat ~ "MUAC category",
    vacc_status = "Vaccination status"
  ),
  intercept = TRUE,
  conf.int = TRUE
)

Characteristic	OR¹	95% CI¹	p-value
(Intercept)	0.01	0.00, 0.04	<0.001
Gender
female	—	—
male	0.89	0.70, 1.13	0.3
Age groups
15+ years	—	—
< 6 months	21.2	4.63, 375	0.003
6 - 8 months	17.5	3.75, 311	0.005
9 - 11 months	19.2	4.06, 343	0.004
1 - 4 years	7.71	1.62, 138	0.046
5 - 14 years	3.31	0.65, 60.5	0.3
Vaccination status
No	—	—
Yes	0.29	0.19, 0.44	<0.001
MUAC category
Green (125+ mm)	—	—
Yellow (115 - 124 mm)	1.39	0.98, 1.93	0.059
Red (<115 mm)	3.82	2.66, 5.45	<0.001
¹ OR = Odds Ratio, CI = Confidence Interval

Time

dat |>
  mutate(
    epiweek = floor_date(date_onset, unit = "week")
  ) |>
  ggplot() +
  geom_bar(
    aes(
      x = epiweek,
      fill = epi_cat
    ),
    position = position_stack()
  ) +
  scale_x_date(
    breaks = "2 weeks",
    date_labels = "%Y -W%W"
  ) +
  scale_fill_manual(
    "Epi status",
    values = c(
      "confirmed" = "#912c2c",
      "probable" = "#c4833d",
      "suspected" = "#edd598"
    )
  ) +
  labs(
    x = "Epiweek",
    y = "N cases",
    title = glue::glue("Epicurve of measle outbreak in Southern Chad"),
    subtitle = glue::glue({
      "{nrow(dat)} cases observed from {min(dat$date_onset, na.rm = TRUE)} to {max(dat$date_onset, na.rm = TRUE)}"
    })
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, size = 5))

By sites

dat |>
  mutate(
    epiweek = floor_date(date_onset, unit = "week")
  ) |>
  ggplot() +
  geom_bar(
    aes(
      x = epiweek,
      fill = site
    ),
    position = position_stack()
  ) +
  scale_x_date(
    breaks = "4 weeks",
    date_labels = "%Y -W%W"
  ) +
  labs(
    x = "Epiweek",
    y = "N cases",
    title = glue::glue("Epicurve of measle outbreak in Southern Chad"),
    subtitle = glue::glue({
      "{nrow(dat)} cases observed from {min(dat$date_onset, na.rm = TRUE)} to {max(dat$date_onset, na.rm = TRUE)}"
    })
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, size = 5)) +
  facet_wrap(~site) +
  gghighlight::gghighlight()

Place

Make a Choropleth map

# clean admin data
dat <- dat |>
  mutate(across(c(sub_prefecture, region), ~ str_to_sentence(.x)))

# count cases by adm2

adm_summ <- dat |> summarise(
  .by = c(region, sub_prefecture),
  n_cases = n(),
  n_deaths = sum(outcome == "dead", na.rm = TRUE),
  cfr = round(digits = 3, n_deaths / sum(outcome %in% c("recovered", "dead", na.rm = TRUE))),
  cfr_lab = scales::percent(cfr, accuracy = .1)
)

# join the count data to the sf

chor_dat <- left_join(
  admin_2,
  adm_summ,
  by = c("adm2_name" = "sub_prefecture")
) |>
  # add AR using population data
  mutate(
    AR = round(digits = 3, n_cases / adm2_pop * 1000),
    label = (paste0(
      "<b>Region:</b> ",
      adm1_name,
      "<br><b>Sub-prefecture:</b> ",
      adm2_name,
      "<br><b>Population:</b> ",
      adm2_pop,
      "<br><b>Attack Rate:</b> ",
      AR,
      "<br><b>CFR (%):</b> ",
      cfr_lab
    ))
  )

leaf_basemap <- function(
    bbox,
    baseGroups = c("Light", "OSM", "OSM HOT"),
    overlayGroups = c("Boundaries"),
    miniMap = TRUE
) {
  lf <- leaflet::leaflet() %>%
    leaflet::fitBounds(bbox[["xmin"]], bbox[["ymin"]], bbox[["xmax"]], bbox[["ymax"]]) %>%
    leaflet::addMapPane(name = "choropleth", zIndex = 310) %>%
    leaflet::addMapPane(name = "place_labels", zIndex = 320) %>%
    leaflet::addMapPane(name = "circles", zIndex = 410) %>%
    leaflet::addMapPane(name = "boundaries", zIndex = 420) %>%
    leaflet::addMapPane(name = "geo_highlight", zIndex = 430) %>%
    leaflet::addProviderTiles("CartoDB.PositronNoLabels", group = "Light") %>%
    leaflet::addProviderTiles(
      "CartoDB.PositronOnlyLabels",
      group = "Light",
      options = leaflet::leafletOptions(pane = "place_labels")
    ) %>%
    leaflet::addProviderTiles("OpenStreetMap", group = "OSM") %>%
    leaflet::addProviderTiles("OpenStreetMap.HOT", group = "OSM HOT") %>%
    leaflet::addScaleBar(
      position = "bottomright",
      options = leaflet::scaleBarOptions(imperial = FALSE)
    ) %>%
    leaflet::addLayersControl(
      baseGroups = baseGroups,
      overlayGroups = overlayGroups,
      position = "topleft"
    )
  
  if (miniMap) {
    lf <- lf %>% leaflet::addMiniMap(toggleDisplay = TRUE, position = "bottomleft")
  }
  
  return(lf)
}

bbox <- st_bbox(filter(admin_2, adm1_name == "Mandoul"))
bins <- c(0, 1, 5, 10, 20, Inf)
pal <- leaflet::colorBin("YlOrRd", domain = chor_dat$AR, bins = bins)
labels <- chor_dat$label |> lapply(htmltools::HTML)

leaflet::leaflet() |>
  leaf_basemap(bbox, miniMap = TRUE) |>
  leaflet::fitBounds(as.character(bbox)[1], as.character(bbox)[2], as.character(bbox)[3], as.character(bbox)[4]) |>
  leaflet::addProviderTiles("CartoDB.Positron", group = "Light") |>
  leaflet::addScaleBar(position = "bottomright", options = leaflet::scaleBarOptions(imperial = FALSE)) |>
  leaflet.extras::addFullscreenControl(position = "topleft") |>
  leaflet.extras::addResetMapButton() |>
  leaflet::addPolygons(
    data = admin_1,
    stroke = TRUE,
    weight = 1.5,
    color = "black",
    fill = FALSE,
    fillOpacity = 0
  ) |>
  leaflet::addPolygons(
    data = chor_dat,
    label = ~labels,
    stroke = TRUE,
    weight = 1.2,
    color = "grey",
    fillColor = ~ pal(AR),
    fillOpacity = 0.3
  )

Fake Measles outbreak analysis

Hugo Soubrier

Epicentre, Paris

last updated 07 Nov 2024

Introduction

Data Cleaning

Age classification

Middle Upper Arm Circumferences (MUAC)

Epidemiological classification

Laboratory data

Person

Demographics

By sites

Age Pyramids

CFR analysis by site

Risks Factor analysis

Time

By sites

Place